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> \ ABSTRACT 

^: 

I We examine emission from a young protostellar object (YPO) with three- 

dimensional ideal MHD simulations and three-dimensional non-local thermody- 
namic equilibrium (non-LTE) line transfer calculations, and show the first results. 
To calculate the emission field, we employed a snapshot result of an MHD sim- 
ulation having young bipolar outflows and a dense protostellar disk (a young 
circumstellar disk) embedded in an inf ailing envelope. Synthesized line emission 
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of two molecular species (CO and SiO) show that subthermally excited SiO lines 
as a high density tracer can provide a better probe of the complex velocity field of 
a YPO, compared to fully thermalized CO lines. In a YPO at the earliest stage 
when the outflows are still embedded in the collapsing envelope, infall, rotation 
and outflow motions have similar speeds. We find that the combined velocity field 
of these components introduces a great complexity in the line emissions through 
varying optical thickness and emissivity, such as asymmetric double-horn profiles. 
We show that the rotation of the outflows, one of the features that characterizes 
an outflow driven by magneto-centrifugal forces, appears clearly in velocity chan- 
nel maps and intensity-weighted mean velocity (first moment of velocity) maps. 
The somewhat irregular morphology of the line emission at this youngest stage 
is dissimilar to a more evolved object such as young Class 0. High angular res- 
olution observation by e.g., the Atacama Large Millimeter/sub millimeter Array 
(ALMA) telescope can reveal these features. Our results demonstrate a powerful 
potential of the synthesized emission of the three-dimensional line transfer to 
probe the velocity field embedded in the envelope, and further analysis will be 
able to determine the precise velocity field to assess the dynamics in the young 
protostellar object to gain a better understanding of star formation. 

Subject headings: stars: formation — ISM : jets and outflows — radiative transfer 
— radio lines : stars - submillimeter 



1. Introduction 



The formation of stars begins with the collapse of a parent dense molecular core. As 
gravitational collapse of a molecular core proceeds, bipolar outflows are expected to form 
in relation to accretion and rotation of the parent core. Observations of the embedded, 
young Class objects revealed that some of them e xhibit infall and rota tional motions in 
emission profiles a s well as bipol ar outfiows (see e.g.. iBelloche et al.ll2002l for IRAM 04191; 
Bourke et allbood for L1041-IRS: Matthews et al.ll20od for Bernerd 1-c). The velocity struc- 
ture in the stage prior to Class objects, before the formation of a protostar, is expected to 
represent the dynamics in the earliest phase of star formation. Therefore, the observational 
studies of the very young stage of star formation, on the way from a prestellar (starless) core 
to a Class object, are quite important in understanding the star formation processes. We 
further can expect a hint of a launching mechanism of molecular outfiows, almost free from 
the propagation effect, in the earliest stage of star formation. 



Before a protostar begins to shine, temperature of the star-forming object without a 
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heating source is as quite low as 10 K. Major observational ways are thus molecular emis- 
sion lines in the millimeter and submillimeter bands. After a protostar forms, the young 
object at the earhest stage of star formation (hereafter "a young protostellar object" or a 
YPO for simplicity) has: 1) a vast amount of envelope matter, 2) complex velocity structure 
consisting of infall, rotation, and outflow motions, and 3) a complicated optical thickness 
structure due to the complex velocity field in the millimeter and submillimeter bands. As 
a result, it is not an easy task to correctly interpret the complex velocity field imprinted 
in the emergent line emission. Furthermore, current observational facilities do not achieve 
sufficient angular resolution to resolve the varying velocity field within the molecular core. 
However, forthcoming high angular resolution observations provided by the Atacama Large 
Millimeter /sub millimeter Array (ALMA) will be able to reveal these internal structures of 
an evolving core, and a significant progress in understanding the star formation process is 
expected. In addition, the launch of the Herschel sateUite will enable observations with 
high frequency resolution in the submillimeter (and shorter wavelength) band. For these 
new observational facilities, detailed modelling and understanding of the radiation field of 
young star forming objects is of crucial importance. In this regard, prestellar cores have 
been studied with one-dirnensional hydrodynamic m odel and radiative transfer calculations 



Pavlyuchenkov et al.ll2008l : iTsamis et al.ll2008l ). Here, we employ a theoretical model 



(e.g., 

constructed with a three-di mensional, high reso lution numerical simulation of the forma- 



tion process of a protostar (IMachida et al.ll2008l ). By combining it with three-dimensional. 



non-local thermodynamic equilibrium (non-LTE) line transfer calculations, we examine a 
characteristic feature in the radiation field of a YPO. In this article, we show the first results 
in formats common to radio observers. Feasibility evaluation of these synthesized emissions 
by the ALMA telescope is another interesting problem, but we do not treat it in detail in 
this article except for simple estimates of the necessary exposure time. 

The outline of this article is as follows. In §2, we describe the models and the set-ups 
of the magnetohydrodynamic and radiative transfer simulations. The results of continuum 
and molecular line emissions are presented in the same formats used in radio observations in 
§3. We make comments with regard to observations with current and future observational 
facilities and discuss them in §4. Summary and conclusions are presented in §5. 



Models and Simulations 



2.1. Magnetohydrodynamic Modelling 



Evolution from a parent dense molecul ar core to a protostar proce eds with a several 



characteristic evolutionary stages (see e.g., iMasunaga fc Inutsuka 



200Cll . and their Figure 
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2b). At the onset of collapse, cooling by dust emission is so effective within a molecular 
core that it contracts nearly isothermally with an initial low temperature. Gravitational 
contraction proceeds and the central density increases, until the core becomes optically 
thick to the dust thermal emission. At this stage the cooling rate is less efficient and the 
collapse becomes adiabatic at the density of n ~ 10^^ cm~^. At this point an adiabatic 
core of a radius i? ~ 1 — 10 AU and a mass ~ O.OIMq surrou nded by the a ccretion shock 



is formed at the center. This core is called "the first core" (ILarsonl Il969l ). The kinetic 
temperature, Tkin, increases as the first core contracts adiabatically. When Tkin exceeds ~ 
2000 K, hydrogen molecules begin to dissociate. Since dissociation is an endothermic process, 
the equation of state becomes soft again, and it induces dynamical contraction of the central 
region at > 10^^ cm~^ (second collapse). Once all the hydrogen molecules are dissociated, 
the equation of state becomes adiabatic again, and eventually a second core (protostar) is 
formed. 



Considering this picture, we employed the MHD modelling of iMachida et al.l (120081) as 
one among those most plausible and based on detailed calculations. iMachida et al.l (120081 ) 
investigated the evolution of a star-forming object from a parent dense molecular core to a 
protostar by three-dimensional resistive MHD simulations. They showed that the magnetic 
force working in an evolving core can naturally drive a slow and wide opening angle flow of 
~ 5 km s~^ and a fast and collimated flow of ~ 30 km s~^ from the first and second cores, 
respectively. In their scenario, these two flows are expected to increase their velocities and 
evolve to a molecular bipolar outflow {v > 10 km s~^) and an optical jet {v > 100 km s^^) 
as observed in more evolved young stellar objects (YSOs). In this article, we follow their 
scenario and simply refer to the low- velocity flow formed in their calculations as an "outflow" . 
Observational study of protostellar objects and outflows in their earliest evolutionary stage 
is indeed quite important for a comprehensive understanding of the star formation process 
as well as for the theoretical modelling. Since the characteristic structures of the youngest 
star-forming objects are small in size, they will be a milestone target of special interest of 
the ALMA telescope. 



Machida et al.l (120081 ) showed that ideal and resistive MHD modelings do not introduce 
a significant difference in the onset of molecular outflows around the first core. In this study. 
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we calculated the cloud collapse based on the ideal MHD equations: 
dp 



+ pyv ■v)v 

dB 



dt 



0, 

-VP - —B X (V X B)- pV0, 

4:71 

V X {v X B), 
4nGp, 



(la) 

(lb) 

(Ic) 
(Id) 



where denotes gravitational potential, and other symbols are used for usual meanings. 

In order t o calculate a model in a wide dynamic range, a nested grid method is adopted 
(see for details iMachida et al.ll2005l ). In this scheme, a rectangular grid of each level has the 
same number of grid points (128^ x 64), and cell width h{l) varies with the grid level I. The 
width of a cell h{l) decreases to a half the width of that of the previous grid level {h{l — 1)). 
A new finer grid is generated whenever the mi nimum local Jeans le ngth Aj decreases below 
8 /i(^max) (thus the so-called Jeans condition of lTruelove et al.lll997l is always satisfied in our 
simulation). We set the maximum and minimum grid level to /max = 10 and /min = 1- The 
grid of / = 10 has a spatial coverage of 250 AU and a cell width of 2 AU. Then our calculation 
can examine the first core and its vicinity with 2 AU at the highest resolution. 

We adopt a critical Bonnor-Ebert sphere of M = 0.6 Mq, R = 2000 AU with a uniform 
kinetic temperature T^in = 10 K as an initial condition. A uniform magnetic field of -B = 
4.2 X 10^^ Gauss threads the initial critical Bonnor-Ebert sphere. Initially the cloud rotates 
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rigidly with dimensionless angular speed tu = ilo/{A-KGpo) ' = 0.3 (where fio = 5.4 x 10 
s~^ is the initial angular speed and po = 3.8 x 10~^^ g cm~^ or uq = 10^ cm~^ is the initial 
central density). We increased the initia l density by 70% to induce gravitational collapse 
(see for more details, iMachida et al.ll2008l ). 



For the thermal evolution, IMachida et al.l (120081 ) employed the barotropic relations P{p) 
ba sed on the equation of s t ate of a one- dimensional radiative hydrodynamic simulation study 
of iMasunaga fc Inutsukal (120001 ). In this paper, we perform a non-LTE line transfer simu- 
lation with a numerical scheme optimized for a uniform Cartesian grid as a postprocessing 
procedure, making use of a snapshot data of the hydrodynamic simulation result of the 
nested grid scheme. In order to maintain a balance between preserving information obtained 
by the nested grid simulation having a fine grid in the innermost region and grasping a 
large scale feature of a whole system, we added a slight modification to the equation of 
state. In particular, we artificially raised the polytropic index after n > 10^^ cm~^ to ^ = 2, 
and stopped further contraction beyond n > 10^^ 



cm ^ (see Eq. [5] of iMachida et al. 



2008 



for comparison). The following set of barotropic relations P{p) was used in the form of a 
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piecewise power-law, 
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where = 190 m s~ is sound speed at T^in 
cm~^) and pd = 3.84 x 10~^° g cm~^ {rid = 1.0 x 10^^ cm~^), respectively. Due to this 
modification, the evolution of low speed outflows from just outside the flrst core can be 
calculat ed for a re l ativel v long timescale. A similar technique was used in the numerical 
study of iTomisakal (120021 ). 

As is described in detail in the next subsection ( §2.2p . we calculated emission from 
the model YPO separately from the magnetohydrodynamic simulation. We made use of a 
part of (/ = 7) snapshot data of magnetohydrodynamic simulation at t = 4200 yrs after 
the formation of the flrst core as an standard input to the radiative transfer code for a 
uniform Cartesian grid. Figure [1] shows the density structure used in the radiative transfer 
simulation. It reveals bipolar outflows extending to about 2000 AU, and a very compact 
and dense flrst core at the center. In the vicinity of the launching region of the outflows, a 
cavity structure exists around the central axis (2;— axis), and a shell of a shape of a letter 
"U" surrounding it. At this stage the major axis of the flrst core is about < 50 AU. In 
addition to the outflows, one can observe a disk-like structure with a radius R < 1000 AU 
an d a scale heigh t < 20 - 300 AU. This disk-like structure also appears in the calculations 
of iMachida et al.l (120081) after a protostar forms. As described above, we did not explicitly 
calculate the evolution of the central region where density exceeds n > 10^^ cm~^. However, 
our long-term calculation effectively enables an examination of further evolution of accretion 
phase to a protostar over a large scale (Figured]), though lack of spatial resolution prevents 
us from observing a protostar formation in a much smaller spatial scale than the minimum 
cell size (~ 10 AU). This disk-like structure is a protostellar disk in its earliest stage, and will 
evolve to a flatter an d rotationally supported circu mstellar disk (hereafter we refer to it as 
a "protostellar disk" ; iHennebelle fc FromangI |2008| ) . To study the formation and evolution 
processes of a rotation-supported circumstellar disk we obviously need calculations of an even 
longer term evolution with higher resolution than the snapshot time adopted here. Since this 
article focuses on the larger scale of young molecular outflows so that emission distribution 
can reasonably be resolved by the ALMA telescope, such an exceptionally high resolution 
simulation is beyond the scope of this article. The protostellar disk is characterized by both 
infall and rotational motions. In order to examine how the emitting and/or absorbing matters 
outside the outflows affect the emergent line intensities, it is convenient to deflne regions of 
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different velocity structure, i.e., by the ratio of infall and rotation velocities. Hereafter in 
this article we denote the least dense gas surrounding the outflows and the protostellar disk 
as an "envelope" , which has only a small ratio of rotational to infall motion and low density 
(n < 2.24 X 10^ cm~'^) in the computational domain. 



2.2. Radiative Transfer Simulations 



We calculate continuum and molecular line emission individually using the snapshot 
data of magnetohydrodynamic simulation described in the previous subsection ( §2.ip . For 
a continuum emi ssion, we calculated thermal e mission from dust grains with opacity of 
Hildebrandl (119831 ) (see also iBeckwith et al.lll99Cll ). For molecular lines, we performed three- 
dimensional non-local ther modynamic equilibrium (non-LTE ) rad iative transfer simulation s 
with an algorithm based on lHogerheijde fc van der TakI (120001 ) (see lWada fc Tomisakall2005l ). 
In the numerical calculation of radiative transfer, rate equations in statistical equilibrium. 



nj ^ Rjj> 

Rjj' 



^jj' + BjjiJ + Cjji J > J', 
Bjj'J + Cjj' J < J', 



(3) 
(4) 



are solved in each cell in the simulation box. In equation (jlj), Ajji and -Bjj' are Einstein's 
coefficients for spontaneous decay and radiative transitions, and Cjji is for collisional tran- 
sition rates from the energy level J to J', respectively. Since almost all the hydrogen is 
in molecular form at this evolutionary stage, we can safely take molecular hydrogen as a 
collisional partner, and then Cjj' becomes 



C 



jj' 



"-Ha 7 J J' 



where is number density of molecular hydrogen, and collisiona l coefficient 'j jji ^ (va) 



is taken from the database tables of Leiden University (LAMDA; iSchoier et al.l 120051 ). In 
order to calculate the mean intensity appearing in equation (jl]) J = (47r)~^ J I^dVL, we 
used a set of sampling rays passing through each cell with a random direction determined 
by Monte Carlo method. Typically the number of sampling rays for each cell is 220 in the 
simulations. As is inferred from Figured! radiation field becomes anisotropic, so we examined 
if the number of ray (A^ray = 220) is sufficient to incorporate the anisotropic radiation field 
by experiments with larger number of A^^ay up to 420. We found only negligibly small 
differences in these results with relative precision less than 0.1%, and the central region 
where coincidence becomes worse than 10% for these runs occupies only a small fraction of 
10~^ for low J lines and a few percent for the maximum J (= 10) in volume. We hence 
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concluded A^ray = 220 is sufficient to obtain reasonably correct results, and below we show 
results of A^ray = 220 runs. 

The specific intensity along the sampling rays are calculated by integrating a standard 
radiative transfer equation, 

(5) 

where is the source function, and J in a cell is obtained from the average of specific 
intensities along the sampling rays passing through the cell. The excitation temperature is 
calculated by iterative calculations of equations (HI) and ([5]) until a self-consistent solution 
of energy level population rij and mean intensity J is obtained. We repeat the iterative 
calculations until the energy level populations rtj conver ge with relative precision of 10 



between the present and former iterations (see for details lYamada. Wada fc Tomisakall2007l ). 
As shown below (§3), optical thickness becomes rather large in most of the lines we calculated. 
To avoid misunderstanding of slow convergence caused by the huge optical thickness as the 
final convergence of self-consistent radiation field, we calculated until the relative difference 
of (i — l)-th and z-th energy level populations fell below 10~^^ or less as an experiment (which 
is taken to be 10~^ as a canonical value in the present paper). Comparison with this longer 
calculation did not find any significant difference from the results with relative precision of 
10^^, and thus we concluded that the level of convergence 10~^ is precise enough to avoid 
false convergence. 

Kinetic temperature of the snapshot result of hydrodynamic simulation is almost isother- 
mal with a low temperature Tkin = lOK, but in a small region at the center Tkm reaches as 
high as 177 K. Then in orde r to include the possible energy lev el cascade from high- J levels 



see for details. Appendix of lYamada. Wada fc Tomisakall2007l ). we examined two choices of 



maximum energy level Jmax = 10 and 16. We found these two runs gave negligible difference 
of 0.1% or less for the levels < J < 10. Therefore in order to save computational time, 
we chose the maximum energy level Jmax = 10 in all the runs in this article. Furthermore 
we confirm our solution by running two cases for each p arameter set with differe nt initial 



energy level populations of LTE and optically thin limit (Ivan der Tak et al.ll2005l ). We ac- 
cept only solutions that meet the criteria that runs of two initial level population coincide 
with relative precision of < 10^^ for both of nj and source function 5*,^. The coincidence of 
S„ is n ecessary because of the rad iative transition terms in equation (jl]), or photon trapping 



effect (lYamada fc Tomisakall2009l ) 



We assume a Gaussian profile (l)y for absorption coefficient only with pure thermal 
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broadening At'th = v^j c^f{^k^L~] 



m] 



(l)[y)dv 



1 



exp 



(6a) 
(6b) 



where z/q is the frequency of local line center, and fc^ is the Boltzmann constant, respectively. 
Doppler effect due to bulk velocity of fluid elements is incorporated through the frequency 
shift of locally estimated vq{v) in each fluid element. Observations of young YSOs such as 
Class objects and molecular cores suggested rela tively small turbulent m o tion compared 



with thermal width froni the line profile analysis (jPi Francesco et al.l 1200 ll : iBelloche et al. 



2OO2I : iBourke et al.l 1200,^ . It is a: 



so observationally suggested that turbulent motion i s 



subsonic in isolated starless cores (IMyersI Il983l : iGoodman et al.l Il998l : ICaselli et al.ll2002l ). 



The snapshot data we use in our radiative transfer simulation is at the evolutionary stage 
on the way from a quiescent starless core to a Class object. Then we neglect turbulent 
broadening in the absorption profile 0(z/). 

The snapshot data from the magnetohydrodynamic simulation cover only a part of the 
computational domain spatially as well as temporally. The magnetohydrodynamic simula- 
tion described in §2.11 is calculated with a fine resolution [dx = 2 AU) by nested grids from 
^min = 1 up to /max = 10, but wc used a sct of / = 7 grid data as a standard. The spatial size 
of the numerical box of / = 7 grid is ~ (2000AU)^x lOOOAU, which covers a whole outflow 
(see Figured]). Precise calculation of emergent emissions from the YPO under consideration 
requires absorption and/or emission of the region outside of the / = 7 grid. However in this 
article, we assume the Cosmic Microwave Background (CMB) radiation of = 2.73 K as 
a background radiation for simplicity, instead of approximate inclusion of outer part. Most 
of the lines in our calculations are optically thick, and thus the effect of this artificial as- 
sumption about background radiation are limited in a very narrow region close to the outer 
boundary of / = 7 grid. In order to estimate the effects of the outer envelope to the radiation 
field of / = 7 grid, we calculated the identical radiative transfer simulation for / = 5 grid 
data which covers a larger volume but has 4 times coarser spatial resolution than / = 7 grid. 
We confirmed that the mean line profiles of / = 5 calculations averaged over the field of view 
do not show qualitative differences from I = 7 calculations especially for mid- and high-J 
lines of SiO molecule. Although a full nested grid radiative transfer calculation is necessary 
to a precise evaluation, we conclude our results of / = 7 grid data do not suffer from the 
significant influence of the outer boundary for SiO lines. For ^^CO and its isotopologue lines 
they turned out to suffer from the outer boundary effect. Thus we will show / = 5 grid data 
results as well as the canonical / = 7 grid data results for ^^CO and its isotopologues. 
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We calculated rotational transitions of the representative molecules used in observations 
of molecular outflows (^^CO, ^^CO, C^^^O and SiO). Among these molecules, we show in 
this article the results of ^^CO (and its isotopologues) as the most popular outflow tracer, 
and SiO as a representative of high density tracer molecule since the average density in 
our snapshot data is high (~ 10^ cm~^). The purpose of this article is to illustrate the 
characteristic features in the emission in relation to the dynamical evolution of YPOs and 
their velocity fields, rather than chemical evolution. So we assume a spatially uniform 
molecular abundance distribution. SiO molecules are in general considered to be formed by 
gas-phase reactions of Si atom desorbed f rom dust grain ma ntles in a shocked gas with high 
velocity > 25 km s~^ (ICaselli et al.lll997l : ISchilke et al.lll997l ). and thus they are regarded as 
a high-velocity shock tracer. However, in the early stage before frozen-out onto dust grains, 
SiO molecule can be formed in gas-phase reactions more than 10 % of the initial Si abundanc e 
in a short timescale 10^ yrs before frozen-out dust onto grains (iLanger fc Glassgoldlll990l ). 
Additionally, observation showed that in some outflows the location of the brightest SiO 
emission does not coincide with the location of the strongest terminal shock (e.g. HH 211: 
Hirano et al.ll2006r). and lower- velocity components than 25 km s~^ have been detected (e.g.. 



Codella et al.l 1 19991 : iGibb et al.l |2004| ). These facts implies that the formation m echanism 



of SiO molecule in protostellar objects remains still ambiguous (lArce et al.l 120071 ). Then in 
order to prevent possible uncertainties in a specific chemical evolution model, we treat the 
abundance distribution as a free parameter in the numerical experiments. Since CS would 
be less sensitive to shock, and it has similar critical densities and wavelengths as SiO, the 
emission of CS lines resemble that of SiO in our calculations with uniform abundance. We 
will show CS results based on the chemistry models in future. HCO^ will also be examined in 
the future article as a charged molecule. We take molecular abundances relative to hydrogen 
molecule as yco = 3 x 10~^ and ysio = 2 x 10^*^ as fiducial values. 



3. Results 

In this section we present our simulation results analyzed in usual methods in radio 
observation. Spectral energy distribution (SED) of continuum emission, maps of synthe- 
sized molecular line intensities integrated along the line-of-sight velocity, position-velocity 
diagrams, channel maps, and line profile maps are presented. While we mainly consider 
high angular resolution observations with the ALMA telescope, we focus on the qualitative 
features in the synthesized radiation field because of the limitation of our simulation set-ups 
(see discussion in §2). That is, we do not include any response of realistic observation in- 
struments by, for example, beam convolution and so forth. We also neglect atmospheric or 
system noises for the same reason. 
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3.1. Continuum: Thermal Emission from Dust Grains 

Theoretically, continuum observation is easier to achieve a higher signal-to-noise ratio 
in high angular resolution observation in comparison to lines. We calculated the thermal 
emission from dust grains in hope of revealing smaller internal structures compared to line 
emission. 

In the snapshot under consideration, the density is as high as n > 10^ cm~^, and we 
can assume a tight coupling between the gas and dust grains. As the kinetic temperature of 
gas is nearly isothermal (Tkin = 10 K) except for a compact region around the first core at 
the center, w e can take the dust temperature Td = Tk in — 10 K. Then dust thermal emission 



is written as (IHildebrandl Il983l : iBeckwith et al.lll990l ) 



I, = 5,{l-exp(-r,)}, (7) 
Ti, = y fimpnKxdl, (8) 

«A = 0.1 i ^ -] cm'g-\ (9) 

where 5^ is Planck function, fi is the mean molecular weight, rrip is proton mass, dl is a line 
element along the line of sight, respectively. In our calculation we take /3 = 2 as a fiducial 
value. Figure [2] shows distributions of thermal emission of = 10 K dust grains in terms of 
brightness temperature Tb = / {2v'^kB)Iv Five panels correspond to different frequencies, 
V = 150, 220, 350, 650, and 850 GHz falling in the band ranges of ALMA receivers. We 
take a single value of inclination angle of the outflow axis to the plane of the sky 6 =30° (see 
a bottom right illustration for deflnition of 6) in Figure [21 We found that average physical 
quantities such as mean intensity do not strongly depend on 6, except for the morphology 
of intensity contours that follow the column density distribution. 

The region bright in dust emission is mainly the central part of the protostellar disk that 
has high column density up to < 10^^ cm~^. Panels (c), (d), (e) of Figure [2] show a weak 
structure of "X" -shape extending from the rim of a box-shape image of the protostellar 
disk. These are base of the slow speed outflows of a wide opening-angle (see Figured]). 
These structures along the outflows are seen clearer in closer view to the edge-on of the 
protostellar disk (i.e. small 6), and are hidden by the geometrically thick protostellar disk 
at > 60° (close to the pole-on view). Note that since the intensity range shown in Figure 
[2] is restricted for a clear view, the contrast between intensities of the outflow cavity wall 
and the protostellar disk is more emphasized compared with the true values. Feasibility 
of detection of the signatures from the outflow cavity in dust emission depends on several 
factors actually, such as sensitivity of a telescope and gas-to-dust ratio in the outflow, and 
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requires more detailed investigation. Figure [2] shows that the brightness temperature of 
dust continuum emission is not strong (~ 0.02 K for 350 GHz), but a huge bandwidth (two 
polarization channels of 8 GHz) of receivers of the ALMA telescope enables the detection of 
these emission easily: for 2 mK sensitivity, the required observation time amounts to about 
2.6 hours for an object at the distance of 140 pc and -30° declination with 0.3 arcsecond 
resolution. For 850 GHz with 1.0 mK sensitivity, required time will be about 8.8 hours, 
which is possible and reasonable to carry out. 

Figure [2] shows that the emission of higher frequency is brighter, whereas equation 
([9]) also implies an increase in optical thickness with frequency. Finite optical thickness 
moderates the increase in intensity (Eq. [7]) compared to the optically thin regime. In Figure 
12] (c), (d) and (e), we indicate the thick {r^, > 1) regions of ~ 20 — 100 AU with solid black 
contours at the center. One can see that though t^, increases with frequency, the region of 
> 1 , in which correction is needed in conversion from to column density, is confined in 
a very small region at the center. As described in §2.H we have to consider the outer region 
surrounding the / = 7 grid employed in our calculation for a more precise estimation. Table 
[1] summarizes the maximum and the mean optical thicknesses in the field-of-view defined by 
the computational domain {1 = 7 grid). It also shows values for (3 = 1.5 and 1. 

Continuum emission shows symmetric distribution with respect to the x = axis [y — z 
plane) and z = axis {x — y plane), reflecting symmetric density structure (Figured]). This 
symmetry is greatly different from line intensity distribution discussed below (Figure^]). 
Panels (c) and (d) of Figure [2] show very weak emission from the narrow structure along 
the outflow axis. Although the dust emission can have information of the density structures 
other than the protostellar disk, it will be quite difficult to detect these weak emissions 
because of much brighter protostellar disk. In spite of our expectation of high angular 
resolution continuum observations, actual observation of dust emission will work as a measure 
of inclination angle of a protostellar disk and configuration of outfiows, rather than examining 
the tiny launching region of outfiows at the center. This is similar to what is done in current 
observations, indicating that we would not take advantage of the high resolution of the 
ALMA telescope. We plot average SEDs for three values of /5 (2, 1.5, and 1) in Figure [3] 
Since the protostellar disk is geometrically thick, mean intensity averaged over the field-of- 
view is almost independent of the inclination angle 9 as well as azimuthal angle around the 
outfiow axis 0, thus multi-band observations will be able to clearly determine /3. Differences 
in SED by (3 appear stronger in longer wavelength regime where optical thickness effect is 
minimal. 



""^We used the ALMA sensitivity calculator for the evaluation throughout the paper 
(http : //www. eso . org/ sci/f acilities/alma/observing/tools/etc/l. 
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3.2. Line Emission: A Probe of Velocity 

In this subsection we show the results of non-LTE line transfer calculations. In contrast 
to the continuum emission, lines become a key to derive information of the velocity field which 
describes the characteristic features in the youngest evolutionary stage in star formation. 

3.2.1. Excitation Status 

In this subsection we examine excitation temperature. In Figure|l]we plot the excitation 
temperatures (Te^{ij) = —huij/kBilia^gi/gj) — \ia{ni/nj)}~^, where gi is statistical weight 
and huij = Ei — Ej) of SiO in each cell as a function of local density. Figure H] shows that 
low- J transition in millimeter band the energy level population is almost fully thermalized 
(Tex — ^kin — 10 K), while high- J transition (upper level of J, Jupp > 4) in submillimeter 
band departs from LTE in the low density regime {n < 10'' — 10^ cm~^). This is due to 
the increasing critical density for LTE with J {ricnt ~ 10^ x cm~^ for SiO rotational 
transitions). The low density region (n < 2 x 10^ cm^^) roughly corresponds to the envelope 
and a part of the protostellar disk, and the high density region {n > 10^ cm~^) concentrates 
to the first core and a part of the midplane of the protostellar disk (see Figured]). Therefore, 
for example, we can expect that J = 4 — 3 (1/43 = 173.883100 GHz) lines will show non- 
LTE effects in the envelope, and J = 7 — 6 (z/76 = 303.9268092 GHz) will show them in the 
envelope and the protostellar disk. We do not show a similar figure of ^^CO: because of lower 
critical density of ^^CO compared to SiO (^icrit ~ 10^ x cm~^), energy level population of 
^^CO is almost completely thermalized up to the highest transition ( J = 10 — 9) in our line 
transfer calculations. 

3.2.2. Integrated Intensity 

For a brief overview of the synthesized line emission distribution, first we show maps 
of integrated intensity of SiO (7 — 6) in Figure O In Figure 0, we overplot red and blue 
components in color lines on the total integrated intensity in gray scale, where intensities of 
total, "red component", and "blue component" are defined with the line-of-sight projection 
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respectively. Four panels in Figure [5] are for different inclination angles [6 =90°, 60°, 30°, and 
0° for panels (a), (b), (c) and (d), respectively). Distributions of red and blue components 
can be described by superposition of rotation of the protostellar disk and the outflows, infall 
and outflow motions. The rotational motion causes symmetric distribution of red and blue 
components with respect to the a; = axis (except for 9 =90° or pole-on view: see panels 
(c) and (d) of Figure [5]), and infall and outflow motions give rise a symmetric distribution 
of them around the center (i.e., /red is almost the same with Jbiue at the location rotated by 
180° around the center: panels (b) and (c) of Figure [5l). Outermost contours surrounding 
the projected outflows in panels (b) and (c) become circular, because the envelope has nearly 
radial infall velocity. 

Next we show the integrated intensity of ^^CO, ^^CO, and C^^O in Figure [61 As is 
discussed in §2.21 these molecules are more easily excited in the envelope compared with 
SiO, and hence the outer envelope surrounding the / = 7 grid cannot be neglected. Line 
transfer simulations of lower I grid data show that the photospheres of lines of these molecules 
are located just inside the radius of the initial Bonnor-Ebert sphere. We show / = 5 grid data 
results in Figure [6] (^ =60°), which cover the whole initial Bonnor-Ebert sphere and thus the 
photospheres of CO and its isotopologues. In these figures we adopted molecular abundances 
y to be 2 X 10-^ 2 x 10"^ and 2 x 10"^ for ^^CO, ^^CO, and C^^O, respectively. Figure 
shows that because of the large optical thickness (see Figure[7] below) , intensity distributions 
have only weak structure compared with SiO (Figure[5]). This result means that, due to the 
high density of the parent molecular core, even the rare isotopologue lines of CO becomes 
optically thick as well as ^^CO. In more evolved system, difference in the optical thicknesses 
in ^^CO and its isotopologues can be used to trace different components, but they are not 
good probe of emission signatures of the very young embedded YPO like that modelled in 
this article. Because the t imescale of CO freeze-out is short in the high density core (~ 10^ 



yrs for > 10^ cm ^: e.g., iFlower. Pineau des Forets fc Walmsleyl l2005l ) . it may be possible 



that significant amount of gas phase CO is locked on the dust grains in the youngest YPO. If 
extremely effective depletion reduces the molecular abundance in the gas phase by a factor of 
10~^ or less in the entire core, the signature of the embedded YPO woul d be easier to appear , 



but such a large factor of depletion is unlikely (see chemistry model, e.g. jAikawa et al.ll200ll ). 
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If the YPO further evolves, emission from the central protostar will heat up the dust grains 
and desorb CO into gas phase again, and furtherr nore the dissociative backgro und radiation 



might prevent CO molecule from freeze-out (e.g.. lJ0rgensen et al.ll2005l . l2006l ). Taking into 
account these factors, the extremely strong CO depletion is unlikely in a real YPO, and 
our conclusion that CO and isotopologue lines are not useful to examine the kinematics of 
a YPO does not change, though absolute intensity and locations of photospheres would be 
slightly affected. 



3.2.3. Average Line Profiles 

Figure [7] plots line profiles of ^^CO emission averaged over the field-of-view for / = 5 
grid calculations. Three panels show different transitions ( J = 1 — 0, 2 — 1, and 3 — 2) for 
a fixed viewing angle 9 =60°. All three transitions display a strongly saturated line profile 
with very weak wing emission {\vr\ ^ 2 km s~^). The strongly saturated profile is due to 
two facts: 1) in our computational domain the density is higher than critical densities for 
these lines (ricrit ~ 10^ x cm~^) and the energy level population is fully thermalized up to 
the maximum J (=10) in our calculation, and 2) resultant optical thickness becomes huge 
[Ty ~ 1500 — 2500 at the rest frequency; note the optical thickness becomes even larger 
if high resolution larger / data are correctly incorporated) as plotted with dashed lines in 
Figure [3 Considering possible depletion of ^^CO molecules on dust grains in the YPO, we 
performed numerical experiments with lower abundances {yco = 3 x 10~^ and 3 x 10~^) using 
the same snapshot data. However, we could not find significant differences in line profiles 
even with a 100 times smaller molecular abundance. These imply that in our calculation 
of a YPO, energy level population of ^^CO is governed by collisional transitions and by a 
relatively small contribution from radiative transitions (oc Bjji) in spite of large optical 
thickness. The same results apply to isotope molecules like ^^CO and C^'^O. Note that as the 
very young protostellar object in our calculation evolves and forms a less dense and larger 
molecular outflow, optical thickness of ^^CO will significantly be reduced and the saturated 
profiles will disappear. The saturated profile in our results, which is not seen in current 
observations, appears because we calculated the emission from the compact high density 
region at the center of the computational domain. 

In Figure [H] we plot average profiles of SiO lines in a similar way but for the canonical 
/ = 7 data results. Since the optical thickness diminishes by a factor of 100 or more compared 
to ^^CO, the average profiles do not saturate but appear as a double-horn shape. Obviously 
SiO lines are better tracers to probe kinematics in the YPO under consideration compared 
with ^^CO (and its isotopologue) lines, besides chemical evolution in the very early stage that 
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has to be further investigated to conclude whether SiO is the best tracer. The double-horn 
profiles appear in all the maps with different inclination angles 6, including pole-on view with 
various degree of skewness in red and blue and varying depth of the dip at Vr = 0. From 
these arguments and Figures [7] and [HI radiative transitions (oc Bjj/) comparatively contribute 
in determining the energy level population of SiO molecule as coUisional transitions. Like 
in the case of ^^CO, we performed numerical experiments of radiative transfer calculations 
with several lower molecular abundances down to ysio = 2 x 10~^°, and confirmed that both 
intensities and profiles of SiO lines change according to the molecular abundance. Our results 
demonstrate that interpretation of subthermally excited lines along with chemical evolution 
should be done carefully. 



3.2.4- Intensity Weighted Mean Velocity along the Lines of Sight 

To display velocity structure clearer, we calculated the intensity- weighted mean velocity 
along the line of sight (or first order moment), 



TfjVrdVr 

(Vr) ^ , (13) 

TijdVr 



and show the distribution of {vr) of SiO(7 — 6) line in Figure[9l As is seen in Figure[Hl in our 
calculations the optical thickness of SiO is also large, though smaller than ^^CO. Then the 
intensity-weighted mean velocity {vr) roughly traces velocity at r ^ 1 region near the outer 
boundary, and hence it does not depend strongly on J. In Figure [H] three panels with different 
viewing angles {6 =60°, 30°, and 0°) are presented. As is seen in Figure [HI (fr) appears 
antisymmetric around the center in panels of 9 =60° and 30° cases (i.e., {vr) — —{vr) at the 
position rotated by 180° around the center). There appears a quite characteristic pattern 
of peaks of (vr) in "S" shape in these panels. These features originate in the rotation of the 
outflows around the ou tflow axis, which are d r iven by magneto-cent rifugal force in the vicinity 



of the first core (e.g., iTomisakal Il998l . l2000l : iMachida et al.ll2008l ). If the outflow having a 



cavity rotates and inclines with respect to the plane of the sky, the total velocity component 
along the line of sight consists of flowing velocity along the outflow axis and rotation around 
the outflow axis. Opposite sign of rotation velocity in two sides of the outflow axis generates 
the gradients of (vr) across the (projected) outflow axis with opposite signs in red and blue 
lobes. If the outflow axis is exactly in the plane of the sky, thus generated gradients in 
(vr) maps appear with the same sign in the two lobes directing upwardly and downwardly, 
because the protostellar disk and the outflows rotate in the same direction. One can see 
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absence of the "S" -shape in the 9 =0° map (panel (c) of Figure[9]), where the outflow axis 
aligns exactly in the plane of the sky. 

In order to demonstrate the effects of rotation of the outflow on line emissions clearer, 
we performed line transfer simulation with a modified snapshot data. We artificially set 
'^<^(= a/'^^x + '^y) = in the snapshot and put it as an input to the line transfer calculation. 
Figure [To] (b) shows the results of f = run of SiO(4 — 3) line. The velocity gradients across 
the outflow axis disappear in the f = run (panel (b) in Figure fTU]) . while the original result 
of SiO(4 — 3) (v<^ 7^ 0) conserves the similar gradients as in FigureO These velocity gradients 
will be able to observationally confirmed with a spatial resolution of ~ 100—200 AU (required 
resolution depends on the inclination angle) and a velocity resolution of ~ 0.1 km s~^. If a 
YPO of this evolutionary stage is at a distance of 140 pc, 100 AU corresponds to 0.71", and 
above structure could be resolved with current instruments such as the Submillimeter Array. 

In the case of CO, because of the locations of photospheres distant from the embedded 
outflow, {vr) appears more weakly than SiO. Figure [11] shows the similar {vr) maps with 
those in Figure [HI but for Z = 5 data calculations. Three panels in Figure [TT] show J = 2 — 1 
transition of ^^CO (y = 2 x 10"^), ^^CO (y = 2 x 10"^), and C^^O (y = 2 x 10^^) for a viewing 
angle 6 =30°, respectively. One can observe that these maps of ^^CO and its isotopologue 
lines quite resemble each other, and they present the weak velocity gradient at the central 
(2000 AU)^ region, of which pattern is qualitatively similar to those in Figure (9] Smaller 
range of (f^) compared to SiO is either because of greater optical depth or smoother density 
structure in the lower resolution / = 5 grid data compared with the / = 7 data. This figure 
again shows that CO and isotope molecules cannot distinguish different components based 
on the difference in optical thickness because of the high density in the model, and puts 
more emphasis on the disability of CO molecules to probe embedded YPOs. 



3.2.5. Position- Velocity Diagrams 

We made position-velocity diagrams from the same synthesized data cubes. First we 
show position- velocity diagrams for ^^C0(2 — 1) and 9 =30° in Figure [T2] for / = 7 data for a 
clear view (we confirmed that / = 5 results are not significantly different from those of / = 7 
calculations). All three panels in Figure [12] show almost no structure (we will discuss the 
origin of a spiky structure in the panel (a) below). These featureless structures in position- 
velocity slices can be easily inferred from the saturated shape of averaged intensity in Figure 
[7] In other words, velocity structures of the hydrodynamic data vanish in the structureless 
Tex distribution and huge optical thickness even in the non-zero velocity part {\vr\ ^ 2 km 
s~^). Panel (c) of Figure [T21 which is a cut passing through the region about 50 AU offset 



- 18 - 



from the center, shows a similar but weaker spiky structure as in panel (a). The origin 
of these spikes in panels (a) and (c) is the same (see below). The cut of panel (b) passes 
through the center, and because of the large optical thickness, it has only a weak velocity 
gradient. 

On the other hand, SiO lines which have much smaller optical thicknesses present a 
variety of structures in position- velocity diagrams due to the non-LTE energy level population 
(FiguredSl SiO(7 — 6) and =30°). The most characteristic features are: 1) some of these 
diagrams show a "gap" around Vr = Q km s~^ (panels (a), (c) and (d)), 2) a pair of cuts at 
the close position (< 50 AU) along the same direction, such as a pair of (b) and (c), presents 
drastically different structures in corresponding posit ion- velocity diagrams. A pair of (d) 
and (e) is another good example. These features are combination of a variety of densities (or 
emissivities) of the (relatively dense) outflows and the protostellar disk that have different 
bulk velocities. Position- velocity diagrams can be used in combination with velocity channel 
maps ( §3.2.61) and line profiles ( §3.2.71) to interpret velocity field. These variety of structures 
in position-velocity diagrams are obviously useful to reconstruct the original velocity and 
density fields from emission lines, although requires high angular resolution observation. 

Finally we mention the origin of the spikes of Aiv ~ ±5 km s~^ in panel (a) in Figures 
[T2] and [131 These arise from the high velocities in the vicinity of the first core at the center 
and thermal broadening. Kinetic temperature in the adopted snapshot is almost isothermal 
of low Tkin = 10 K, but in the limited region close to and in the first core Tkin rises rapidly 
as high as 177 K at a maximum, which corresponds to the thermal width Af^h ~ 2 km 
s~^. In our numerical model, the magneto-centrifugal force associated with the "hour-glass" 
shaped magnetic fie lds formed around the rotating first core d rives the outflow of a wide 



opening angle (e.g., iTomisakal Il998l . l2000l : iMachida et al.l 120081 ). In this model, magnetic 



field is straight on the axis (the outflow axis) and thus the magneto-centrifugal force 
does not work there. The velocity near the 2;— axis is not outward, but inward. At the 
center where gravitational potential depth is at its maximum, the inflow motion takes its 
maximum velocity as high as 1.8 km s~^. The combination of a wide thermal width and the 
highest velocity (in the whole snapshot data) results in spikes in panel (a) in Figures [12] and 
[13] with a very large velocity range in the position- velocity diagrams. Because of the large 
thermal broadening, velocity range at these spikes becomes larger than the maximum bulk 
hydrodynamic velocity of the input data. To confirm the origin of these spikes, we performed 
line transfer experiments with the same snapshot but its kinetic temperature is set to be 
10 K isothermal. The spikes disappear in these experiments, as expected. 



- 19 - 



3.2.6. Velocity Channel Maps 

As an alternative expression of the velocity field imprinted in line emission, we show 
the velocity channel map of SiO(7 — 6) line in Figure [TH In Figure [TH extended diffuse 
components appearing in the velocity range -0.6 km s^^ ^ ^ 0.6 km s~^ are emission 
from the rotating protostellar disk. The accreting envelope extending further outside does 
not contribute in the channel map, since its density {n ~ 10^ cm~^) is smaller than critical 
density of SiO(7 — 6) {ncnt — 1-3 x 10^ cm~^ at Tkm = 10 K). The density in the snapshot 
has a cavity-like structure and a slightly denser fan-shape structure surrounding it (see 
Figured]). However in the channel map (Figured!]) the fan-shape structure close to the 
outflow launching region is hardly observed because of the geometrically thick protostellar 
disk, even in the edge-on view. Only in a very limited range of velocity channels and viewing 
angles can we see the fan-shape structure. This fan-shape structure is due to the magneto- 
centrifugal force driven flow from the rotating first core, and does not come from entrainment 
of the core matter by a high-velocity jet. Our results mean that close attention should be 
paid in inferring a driving mechanism of outflows from observed channel maps especially 
when the target protostellar object is at the young evolutionary stage. 

One of advantages of velocity channel maps is that it can demonstrate the rotational 
motions of the outflows more clearly compared to other diagrams. For a comparison we 
plot a channel map expected for the identical dataset except for f = in Figure dSI (see 
§3.2.4p . Comparison with Figure dH clearly shows the existence or absence of rotation of 
the outflows. High angular resolution observations with the ALMA telescope will reveal 
these complex structures of rotating outflows coupled with a protostellar disk. Note that 
further imaging simulations are necessary to examine whether the diffuse emission from 
the protostellar disk can be correc tly imaged with ALMA or other interferometers (e.g.. 



Kurono. Morita &: Kamazakil l2009l ). Tentative experiments showed that the total power 
array of the ALMA telescope is crucial to reproduce the diffuse emission in Figures dH and 
[T5] (Kurono & Yamada, private communication). 



3.2.7. Line Profile Distributions 

Finally we present an example of line profile distribution map in Figure [T6] (SiO (7 — 6) 
of 6' = 30° view). Figure [T6l shows that subthermal lines have double-horn or multi-peak 
profiles almost over the entire field-of-view. A nearly flat profile at the map center comes 
from a combination of large bulk velocity near the first core and a broad thermal width 
in the region ( §3.2.5p . Blue component is relatively brighter on average (Figure[H]), but 
opposite trend (brighter red component than blue one) appears locally. This average blue- 
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ward skewness would largely come from the overall contraction of the core (IZhoulll995l : [Evans 
1999h . 



A protostellar object at the very young evolutionary stage has a quite complicated ve- 
locity structure composed of accretion of the parent molecular core, rotation, and outflow 
motio n and so forth. Therefore, as often referred in s tudies of spherical starless or prestellar 
cores (jPavlyuchenkov et al.ll2008l : iTsamis et al.ll2008l ). the relation between emission features 
regarded as signatures of infall motions (e.g., blue-red asymmetry in line profiles) and actual 
infall motions is not necessarily well-established. For example, in our calculations, the snap- 
shot data has an infall motion in almost entire computational domain and optical thickness 
large enough for self-absorption , nevertheless loca l enhancement of red intensity over blue is 
observed in some lines of sight. lAdelson &: Leungl (119881 ) found that rotation and expansion 
motions can show a similar double-horn profile with a brighter red component compared 
to a blue one. Interestingly, our simulation results show double-horn profiles even in the 
outer part, though weak, where the optical thickness is not very large. The asymmetric 
double-horn profile is basically a result of combination of several factors such as: 1) blue- 
enhancement by asymmetr ic self-absorption due to the increasing Tex toward the center in 
infall motion (jEvanslll999l ). and 2) bipolar velocity (infall, outflow, and rotation) reflected 
in the optically thin regime. An additional numerical experiment is performed in order 
for a simple examination of emissions from different velocity components. We constructed 
two modified input data that have 1) radial infall velocity field outside the outflow regions 
(v = Vj. = mm{v • Bj., 0)ej., where e,. is a unit vector in the radial direction; case A), and 2) 
the outflow region velocity only {v = (i^out • ^z)^z'i case B). In case B, the outflow region 
velocity Vout is defined as the velocity having If^l > 0.5 km s^^: this criterion turned out to 
be able to effectively capture the outflow region after several experiments. Figure [17] shows 
the results of line transfer simulation using these artificially constructed data. Panel (a) of 
Figure [T7] shows the result of case A. It shows blue-red infall asymmetry almost over the 
entire field-of-view and confirms this asymmetry in Figures M and [TH] originates in overall 
infall motion of the core. Panel (b) of Figure [T71 on the other hand, represents almost sym- 
metric double-horn profiles especially in the vicinity of the projected outflow axis. Since in 
this result the infall velocity is artificially fixed to be except for that on the z— axis (see 
§3.2.5p . this symmetric profile would come from the bimodal outflow motion itself, and the 
skewed profiles close to the z = plane would be due to the protostellar disk. In the full 
data set simulation, these components couple nonlinearly with structured emissivity, optical 
thickness, and excitation temperature in different velocity and position. Further detailed 
analysis is necessary to disentangle them in a qualitative way. 
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3.3. Velocity Structure and Line Emission 

As displayed in the previous subsection §3.21 line emissions from YPOs are quite compli- 
cated. In Figure [18] we plot relations between bulk velocity and local density of the adopted 
snapshot data in the two-dimensional probability distribution function (PDF). Figure [T8] 
shows that: 1) distribution of and \v\ with respect to density is qualitatively similar, 
and 2) either of f = or \v\ = components are not significant in mass. The velocity of 
the magneto-centrifuga l force driven outflows are roughly equal to the Kepler speed (e.g., 



Kudoh fc Shibatalll997l ). and then in the adopted snapshot, Vj-ot, "Winf and Vout (rotation, infall 
and outflow speeds) become comparable. Hence the first feature in Figure [TSl similarity in 
the distributions of and represents comparable speed of v^ot, finf and Vout- Taking into 
account the symmetric structure of the input snapshot with respect to the z = plane and 
the geometrically thick configuration of the protostellar disk, it explains a complex velocity 
structure that multiple velocity components with similar norms but different signs easily ap- 
pear in the field-of-view (or on a line-of- sight) whichever inclination angle aligns the object. 
In a more evolved object such as seen in current observations, the velocity structure becomes 
simpler mainly because a protostellar disk will evolve to a thinner rotationally supported 
circumstellar disk with relatively small spatial extent compared to the outflows. Hence a 
complex velocity feature is a manifestation of dynamical evolution of the star-forming object 
in the young evolutionary stage, when the outflows are embedded in the accreting envelope 
and the thick protostellar disk. Therefore such complexity in emission induced by complex 
velocity structure is a characteristic of the early stage of star formation. 

We plot the relation between bulk velocity of the snapshot data and excitation temper- 
ature obtained in line transfer calculations in Figure [T^ for SiO(2 — 1) and (7 — 6) lines as 
two-dimensional PDFs. Figure [12] shows that Tex has a wide scatter as a function of veloc- 
ity, and hence self- absorption will be possible at various velocities. Particularly, a relatively 
large scatter of T^x around f ^ = and Vz = suggests a large possibility of self-absorption 
at zero velocity as well. Either a small amount of matter with v = (Figure [T8|) and high 
probability of self-absorption at f ^ = or f ^ = (Figure [T^ can weaken the line emission at 
Vr = compared to the other velocity components. This causes the gap at f ^ = in position- 
velocity diagrams in Fig ure [T3l The gap structur e in position- velocity diagram were al ready 



seen in observation (e.g. iHogerheijde et al.lll998l : [Williams et al.ll2003l : [Lee et al.l [20061 ). On 



the other hand, energy level population of ^^CO is almost fully thermalized in the nearly 
isothermal snapshot, thus T^^ — Tkm — lOK and ^^CO lines do not show a gap in position- 
velocity diagram (Figure [T2]) . Since neither Figures [T8] nor [T9] include spatial correlation, 
they cannot explain all the features of emitted lines. All the arguments above are statistical. 
Further detailed analysis of spatial structures will appear in the following article. 
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Discussion 



The early phase of protostellar objects and their evolution in its beginning are quite 
interesting and importa nt problems in the study of s tar formation. For the formation of low- 
mass stars, a model of lAdams. Lada fc Shul (119871 ) of spectrally categorized young stellar 
objects from Class 1 to Class 3 has been widely accepted. While multi-dimensional studies 
using numerical meth ods have constructed r ealist i c evolutionary models of molecular cores 



to protostars fsee e.g.. lTomisakalll998l. l2000l . l2002l : iBanerjee fc Pudritzll2006l : iMachida et al. 
2006 . 2008 : Hennebelle fc Fromangj 2008 ). observational studies have proceeded on the basis 



of classical spherical models. One of the difficulties in observational studies of the early stage 
of star formation is its deeply embedded signature in the parent molecular core. Observa- 
tions have found a riumbe rs of non-radial velocity fields in mostly spherical molecular cores 
(e.g., iBelloche et al.ll2002l ). but they have been interpreted in terms of a simple superposi- 



tion of arbitrarily introduced rotation or outfiow/infiow motions within the framework of 
classical spherical models. High angular resolution observation facilities like the ALMA tele- 
scope are expected to provide information of velocity fields even in small structures within 
a protostellar object. However, as shown in the previous sections, the signals of emission 
from a YPO deeply embedded in the envelope are difficult to correctly decipher even with 
the high angular resolution of the ALMA telescope, because of the complex velocity field 
and resultant complicated optical thickness structure. Furthermore identification of outfiow 
itself would be far more difficult compared to the case of more evolved outfiows currently 
observed (FigurefT^. Therefore, careful investigation of how the true structures of objects 
are refiected in emission from them is necessary. 

Recently the Spitzer Space Telescope found numbers of very low luminosity objects 
(VeLLOs) in the mid-infrared band. The distribution of molecular gas accompanying these 



objects and the existence or absence of outfiow are quite uncertain (IBourke et al.l 12006 



Yamada & Ohashi, private communication). If these VeLLOs are indeed related to the 
very early phase of star formation, irregular morphology of molecular emission might also 
be a natural consequence (e.g., FigureEj). If future high-angular resolution observation 
of molecular emission of VeLLOs indicates an outfiow signature, comparison of our results 
would provide a valuable insight on the earliest phase of star formation. Among VeLLOs, 
for example, a molecular core L1521F in Taurus has a high central density, infall asymmetry, 
and is chemically evolved. Since the discovery of a weak infrared source by the Spitzer 
telescope, it is considered to be a n object on the evo lutionary way from an evolved starless 
core and Class protostar stage (IBourke et al.ll2006l ). thus it would be a good candidate to 
be studied in the future high-angular resolution observation with the ALMA telescope to 
test the prediction of our paper. A more evolved VeLLO, IRAM 04191+1522 also in Taurus 
will be another good target of the ALMA telescope in the study of the early stage of star 
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formation. Besides VeLLOs, it might be possible that rotation signatures hke Figure [H] were 
considered to be a manifest of oth er kinds of velocities in so me of the observation of irregular 
outflows in young YSOs (see e.g., iDi Francesco et al.ll200ll ). 



4.1. Relevant Radiative Transfer Simulation Studies 



Recently line transfer simulation studies of star formation have increased in number. 
Many works on the early phase of star formation focus on c hemical evolu t ion in their one- 
dimensional modelling rather tharion dynam ical evolution. iTsamis et al.l (120081 ) employed 
the inside-out collapse solution of IShul (119771 ) and calculated chemical evolution along with 
line transfer in collapsing cores. They constructed a virtual 45-m class telescope and per- 
formed pseudo-observation of their simulation results. iPavlyuchenkov et al.l (120081 ) system- 
atically examined the effects of rotation and infall motions on emergent line emission with 
their one- dimensional modelling. Complex velocity fields (infall, rotation and outflow) in the 
early protostellar objects affect the emission line profiles in a quite non-trivial way ( §3.2.7p . 
Generally blue-skewed asymmetry and inverse P-Cygni profiles are regarded to evidence an 
infall motion as well-understood mechanisms. However, quantitative understandings of the 
effects of gas motion on line emission is not significantly well developed to solve an inverse 
problem. In other words, we currently still cannot reconstruct the velocity field from the ob- 
tained line profiles to a satisfactory degree for a comparison with theoretical models, even for 
a spherical gas. These works are experimental rather than applicable for asymmetric objects 
in real space, nevertheless they will form a powerful set of building blocks in interpreting 
more realistic simulation results like ours and observational results. 

Many numerical works have been done in order to derive physical quantities of outflows 
i n obs ervation data so that models can fit the observational results. iMevers-Rice fc Lada 
(11991 ) performed numerical experiments of a simple concave model outflow . andlRawlings et al 



(2004) succeeded in constructing a best fit structure to the results of iHogerheijde et al. 



(119981 ) . that consists of outflow cones and a surrounding envelope with their two-dimensional 



radiative transfer calculations. These studies employed minimal models to reproduce the 
temporal structure imprinted in observations. On the other hand, in more realistic calcula- 
tions presented in this article, the most difficult factor in the identification of the emission 
components from the outflows is the contamination from the geometrically and optically 
thick protostellar disk (Figure fT^ . In a more evolved system, outflows cover wider spatial 
extents than the protostellar disk, and the protostellar disk will evolve to a flatter circumstel- 
lar disk, so that the outflows are observable without shielding of a thick disk. In other words, 
the youngest outflows would not appear similar to those seen in current observation. Inter- 
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ferometric observational studies proposed a disk/outflow model based on multi-components 
appearing in their data, but with an assumed geometrically thin disk. Possible effects on 
the emission by a thick protostellar disk on the evolution toward a circumstellar disk is first 
mentioned in th i s article. A thick protostellar disk appears in a number of numeri cal studies 
(iTomisakal llQQsl . boool . I2OO2I : ISaneriee fc Pudritzl bood : iMachida et aP bood . boosh . Detailed 
study of emission is necessary for a verification of these results in future observations. 



4.2. Launch of Outflows 



Our calculations have a relation to outflow launching mechanisms. One of the charac- 
teristic features of outflows driven by magneto-centrifugal force is its rotation around the 
flow axis. In contrast, if the outflow is accelerated by the momentum transfer from the jet 
(so-called "entrainment" mechanism), the molecular outflow would rotate very slowly. This 
is expected by the fact that, in the entrainment mechanism, the angular momentum is not 
effectively redistributed from the jet which has a small angular momentum to the molecular 
outflow which has a large moment of inertia. Our line transfer calculations have shown that 
the rotation o f the outflows appears in velocity channel maps, position-velocity diagrams and 
maps of {vr). iLaunhardt et al.l (120091 ) recently found a gradient of {vj.) supposed to present 
rotation of an outflow form a transition object from Class 1 to Class 2 (or younger), CB26 
in their IRAM PdBI observation. In addition to an evolved system like CB26, if rotation of 
outflows in younger objects will be discovered, they will be a good constraint for the models 
of outflow launching mechanisms. The younger the object, the tighter the constraint will be, 
though a younger system is embedded in an envelope of parent molecular core matter, and 
its definite detection is more difficult. We found a signature from rotation of the outflows in 
velocity channel map, position- velocity diagrams and distribution of (f^). Further analysis 
of the line profile distribution as well as those in Figure [17] will provide clearer indications 
of the rotation signatures, and will appear in a future article. 



In the evolution model we employed in this article (jTomisakalll998l . l2000l . l2002l : lMachida et al 



20081 ). outflows are driven just outside a first core. The strongest support for their model 
will then come from the emission from the first core. However, we cannot derive any strong 
conclusion about this, since the numerical modeling adopted in this article cannot resolve a 
geometrically thin shocked region at the surface of the fi rst core, and canno t follow the shock 
chemistry inside such a thin region. lOmukail (120071 ) and lSaigo fc Tomisakal (120091 ) calculated 
emission from a thin cooling shock surrounding a first core. They consistently examined 
emission from the radiative shock and absorption / re-emis s ion in the surrounding envelope 
by sacrificing the detailed dynamical evolution. lOmukail (120071 ) found that reprocessings 
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in the thick envelope is significant in a spherical cloud, whereas ISaigo fc Tomisakal (120091 ) 
showed that rotation flattens the cloud and reduces the column density along the rotation 
axis by more than one order of magnitude. Our three-dimensional calculations of longer 
wavelength emission, mainly for the ALMA telescope, showed that line emission can be seen 
thro ugh the wider v eloci ty range due to the non- spherical motions and geometry. Different 
from lOmukail (120071 ) and lSaigo fc Tomisakal (120091 ) . our calculations do not explicitly include 
an outer absorbing envelope ( §2.2p . Though we conclude that the outer envelope would not 
qualitatively alter our results of SiO lines of / = 7 grid by a simple survey of line profiles in 
/ = 5 calculation, these effects in outer envelope absorption and optical thickness should be 
included in an advanced numerical scheme in the future. 



4.3. Comments on the Real Observations 

Finally we briefly comment on the real observations especially with the ALMA telescope. 
Our calculations are the appropriate means for comparing the results of MHD calculations 
with observations since we calculate the radiation explicitly. This article focuses on the study 
of characteristic features of radiation field of a YPO, and does not include ar iy limitations 



from observational instruments such as finite beam size (compare with e.g., iTsamis et al. 



20081 ). Although our results are ideal, we can make a simple estimate of the necessary 
observation time by the ALMA telescope using the synthesized emission presented in this 
article. If we assume a YPO at the distance of 140 pc and declination S = —30°, typical size 
of the characteristic features in the synthesized maps (Figures |5l [6], [14] and [15]) is > 10 AU, 
or ~ 0.1" in angular size. From FigurelH the brightness temperature of SiO(7 — 6) line is ~ 3 
K, and line width is ~ 2 km s^^ . Then if we invoke the velocity resolution to be 0.1 km s~^, 
which is used in our calculation, and 0.3 K sensitivity, the necessary exposure time amounts 
to ~ 5100 seconds or ~ 14 hours. For the continuum, the required exposure time reduces 
to ~ 2.6 hours for the same object but 2 mK sensitivity because of a large bandwidth of 
the receivers installed on the ALMA telescope ( §3.ip . It is noted that one should be careful 
whether diffuse emission from the protostellar disk in Figures [TJ] and [TB] can be correctly 
reproduced in interferometer observations or not. The tentative results of such imaging 
simulation for the ALMA telescope suggest that in order to reproduce the diffuse component, 
the total power array is crucially important, and the above estimates assume the use of the 
full array (12 m, 7 m, and the total power array; Kurono & Yamada, private communication). 
Lower J lines, of which critical densities are lower than the average density in the protostellar 
disk and the envelope, appear even extended than Figure [H] Besides the limitation that 
comes from radiative transfer scheme ( §2.2|) . imaging simulation for i nterferometers will be 



necessary for direct comparison of our results with actual observations (iTakakuwa et al.|]2008 
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Kurono et al. 2009) 



5. Summary 

In this article we surveyed the characteristic features of the radiation field from a dy- 
namically evolving object at an early stage of star formation, on the way from a prestellar 
(starless) core to a Class object. We have presented the most realistic calculation results 
obtained to date from a combination of three-dimensional MHD simulation and non-LTE 
line transfer calculations. We calculated the dust thermal emission and rotational lines of 
^^CO (and its isotopic molecules) and SiO molecules in the snapshot data consisting of the 
compact outflows {R ~ 1000 AU), the protostellar disk and the envelope. We found that: 

1) Dust thermal emission traces a dense mid-plane in the protostellar disk rather than 
very small structures in the launching region of the outflows. A very weak fan-shape signature 
along the base of the outflows appears in viewing angles close to edge-on, but it is likely to 
be overwhelmed by the much brighter protostellar disk. 

2) Because of the high density {n > 10^ cm~^) in our snapshot data, ^^CO and its iso- 
topologue lines are fully thermalized and are not suitable to probe velocity field. Our results 
are in apparent disagreement with current observation of more evolved outflows, but can 
be reasonably understood in terms of the different densities at different evolutionary stages. 
This implies that the energy level population of CO is governed by collisional transitions, 
and the molecular abundance does not significantly affect the results. We experimented with 
lower molecular abundances to mimic possible depletion of ^^CO molecules, and found no 
significant differences in the final results. 

3) On the other hand, SiO lines are not completely thermalized because of their higher 
critical densities [ricnt ~ 10^ x cm~^), and are more suitable for probing the velocity 
structure in contrast to CO. Line profiles of SiO are characterized by a double- horn shape over 
almost the entire field-of-view, and position- velocity diagrams showed a variety of structures 
strongly dependent on locations of different cuts. These structures arise from a complex 
distribution of optical thickness and non-LTE energy level populations induced by several 
velocity components, such as outflow, inflow, and rotation. In a dynamically evolving system, 
these velocity components have similar speeds but different directions, thus any line-of-sight 
can pick up multiple velocity components of similar norms but different signs irrespective 
of viewing angles. In order for a correct interpretation of such complex velocities from line 
emission, further detailed analysis is necessary. 



4) In the adopted model in this article, the outflows are driven by a magneto-centrifugal 
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force (iMachida et al.l l2006l . |2008| ). One of the characteristics in magneto-centrifugal force 
driven flows is their rotation around the axis. Its signature can appear in intensity-weighted 
velocity maps and velocity channel maps of the synthesized radiation, and should be ob- 
servable with the ALMA telescope with a reasonable time (~ 14 hrs for SiO(7 — 6) line and 
~ 2.6 hrs for continuum at 350 GHz). Se veral mode l s hav e been proposed for launching 



mechanism of molecular outflows (see e.g., lArce et al.l 120071 ). but are still in dispute. Our 



results showed a clue to identify the signature of magneto-centrifugal force driven models 
(rotation of the outflows: Figured!]), and will be an important step towards understanding 
of the early stage of star formation. 

Based on the limitation of the outer boundary for the radiative transfer calculations, 
our results should be taken as rather qualitative, and not quantitative for comparison to real 
observations. Further studies will be needed for a confronting theoretical model prediction 
for the ALMA telescope and for revealing, in particular, the connection between velocity 
and emission lines. 

We thank T. Matsumoto and T. Hanawa for contribution to the nested grid code used in 
hydrodynamic simulation. M.Y. thanks for intriguing and inspiring discussion and comments 
of N. Hirano, C.-F. Lee, N. Ohashi and M. Momose, and R. Taam for his critical reading 
of this article. Numerical computations were partly carried out on VPP5000 and Cray XT4 
at the Center for Computational Astrophysics (CfCA), National Astronomical Observatory 
of Japan, and this research was supported in part by Grants-in-Aid by the Ministry of 
Education, Science, and Culture of Japan (17340059, 16204012,18740104). 
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Table 1: Mean optical thickness of dust emission averaged over a field-of-view {ti) and the 
maximum optical thickness in the field-of-view r^^^max for ^ = 30° view. Values in this table 
do not strongly depend on (see Figure[3]). Three cases with different power- law indices of 
opacity /? = 2, 1.5, and 1 (Eq. [9]) are listed. 
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Fig. 1. — Density structure of the snapshot data used in radiative transfer simulation. The 
density at x — z plane at ?/ = is displayed. Bipolar outflows extend along the 2;-axis 
{R < 1000 AU), and the first core is seen as a quite compact dense region {R < lOOAU) 
at the center. Additionally a geometrically-thick protostellar disk appears as an extended 
region with a scale height ^ 200 — 300 AU. In the vicinity of the launching regions of the 
outflows, the density structure shows a cavity surrounding the central axis along the 2;— axis, 
and a shell-like dense regions of the shape of a letter "U" surrounding them. 
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(a) 150GHz 




Fig. 2. — Maps of thermal emission from = lOK dust grains. Five panels correspond to 
different wavelengths in the millimeter and submillimeter bands. All these maps have an 
inclination angle 6 =30° (see the bottom right illustration for definition of the inclination 
angle 9). Black solid lines at the center (i? ~ 20 — 100 AU) of the panels (c), (d) and (e) 
indicate regions of which optical thickness exceeds 1. Correction of finite optical thickness 
for conversion from Tf, to column density seems unimportant except for an innermost region 
even in the highest frequency in the submillimeter band in our simulation, but thick regions 
broaden their areas if the outer part of the simulation domain is included. 
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Fig. 3. — SED averaged over the field-of-view defined in the computation domain. Blue hues 
are for (3 = 2, red hues are (5 = 1.5, and yellow lines are /? = 1, respectively. For a value of 
/3, solid line indicates the viewing angle of 6 =90° (pole-on view), dotted line is for 6 =60°, 
dashed hne is for 9 =30°, and a dot-dashed line is for 9 =0° (edge-on view). 
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Fig. 4. — Excitation temperature of SiO line in each cell is plotted as a function of density 
of the cell. Thin horizontal lines at Tex = 10 K indicate a reference for thermalization in an 
isothermal gas of Tkin of lOK. Low excitation transitions in millimeter bands {J = 1 — 0, 2 — 1) 
are almost fully thermalized, but the excitation conditions of high excitation transitions of 
which critical density for LTE rapidly increases with J (ncrit ~ 10^ x cm"'^), depart from 
LTE in the low density regime (n < 10^ — 10^ cm~^). 
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Fig. 5. — Maps of integrated intensity of SiO(7 — 6) line (gray scale). Four panels correspond 
to different viewing angles: panel (a) is 6 =90°, (b) is 6 =60°, (c) is 9 =30°, and (d) is 6 =0°. 
Red and blue contours indicates integrated intensities in f ^ > (red) and f ^ < (blue) 
components, respectively, with equal contour spacing of 1.05 K km s^^. 
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Fig. 6. — Maps of integrated intensity of ^^CO, ^^CO, and C^^O (1 — 0) lines taken from the 
/ = 5 data calculations. Because of the use of / = 5 grid data, the area of the field-of-view is 
16 times larger than that in Figure [51 Contours are drawn with equal spacing of 1.05 K km 
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Fig. 7. — Line profiles of ^^CO averaged over the field-of-view (solid lines) and average Ti, 
profile (dashed lines) as a function of line-of-sight velocity for 6 =60° view for I = 5 grid. 
Panels (a), (b), and (c) are for different transitions, J = 1 — 0, 2 — 1, and 3 — 2, respectively. 
Because of the huge optical thickness, line profiles take a strongly saturated profile with very 
weak wing signatures {\vr\ ^ 2 km s~"^). 
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Fig. 8. — The same figures as Figure[71 but for SiO. Panels (a), (b), and (c) denote J = 2 — 1, 
4 — 3, and 7 — 6 transitions, respectively. Mean optical thickness reduces by a factor of ~ 100 
compared to ^^CO, and line profiles take a double-horn shape skewed to the blue side. 
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Fig. 9. — The velocity first moment ((fr)) maps of SiO(7 — 6) line. Panels (a), (b), and (c) 
are different inclination angles {6 =60°, 30°, and 0°, respectively). 
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Fig. 10. — The velocity first moment (vr) maps of SiO(4 — 3) line for 6 =30° view. Panel (a) 
is taken from the original result with 7^ 0, and panel (b) is a result of f ^ = experiments. 
Symmetric and anti-symmetric distributions of peaks are due to rotation of the outflows and 
the protostellar disk. Similar trend is observed in all the lines we calculated. 
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Fig. 11. — The velocity first moment ((fr)) maps of ^^C0(3 — 2) and its isotopologue line 
for / = 5 calculations. Panels (a), (b), and (c) are ^200(3 - 2), ^^C0{3 - 2), and C^^O, 
respectively, and have the same viewing angle 6 =60°. The global symmetric gradient of (f,.) 
would originate in the rotation of the parent core. 
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Fig. 12. — Position-velocity diagrams of ^^C0(2 — 1) line along three cuts indicated on the 
top- left panel ((fr)) for 9 =30° view (data are taken from / = 7 grid calculation for a clear 
view). Width of position- velocity diagrams are scaled to the length of the cuts from the 
center. Solid lines indicate {vr) along each cut. 
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Fig. 13. — Position- velocity diagrams of SiO(7 — 6) line. Horizontal width of diagrams are 
similarly scaled with the cut length on the {vr) map at the top- left. Compared with Figure 
fT2l SiO line shows a variety of structures in these position-velocity diagrams due to the 
non-LTE population. 
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Fig. 14. — Velocity channel map of SiO(7 — 6) line for 9 =30° view. The relative velocity 
in unit of km from the rest of the object is labeled at the top-left of each map. Diffuse 
components appearing in \vr\ < 0.6 km s~^ are emission from the protostellar disk. Rotation 
of the outflows and the protostellar disk is clearly observed. 
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Fig. 15. — Same as Figure UM, but for a f <^ = run. No rotation signature from the outflows 
and the protostellar disk is observed. 
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Fig. 16. — Line profile distribution of SiO(7 — 6) for a 6 =30° view are sampled by ~ 300 AU 
spacing over the entire field-of-view. Solid lines denote intensity profiles, and dashed lines 
are optical thickness on the same line of sight. Line profiles in double-horn shape are almost 
ubiquitously seen in this view with varying relative ratio of red and blue peaks. 
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Fig. 17. — The same distribution maps of line and optical thickness profiles as Figure | 
but for the experimental results for case A (panel (a): radial velocity only) and B (panel (b): 
z-component of Vout only). See text for details. Panel (a) shows the blue-skewed double-horn 
profiles almost in the entire field-of-view, and panel (b) represents symmetric double-horn 
profiles in the outflow region. 
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Fig. [T7] cont. — Line and optical thickness profile distribution for the case B experimenet 
(see text). 
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Fig. 18. — Velocity structure of the adopted snapshot data in two-dimensional PDF. Panel 
(a) shows the relation between and density, and panel (b) shows that of and density 
in both color scale and contours (0, 10°, 10^, 10^, 10^, and 10^). 




Fig. 19. — Relations between bulk velocity structure of the snapshot and excitation temper- 
atures of SiO(2 — 1) and (7 — 6) lines are shown in two-dimensional PDFs (solid contours 
are 0.0, 0.75, 1.5, 2.25, 3.0, 3.75). Panels (a) and (b) describe the relation of Tex and v^, 
and panels (c) and (d) describe that of Tgx and Vz- Higher J line shows a broader scatter for 
either of v^f, and Vz-, and has a higher possibility of self-absorption compared to lower J line. 
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